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Recently, crust cooling times have been measured for neutron stars after extended outbursts. 
These observations are very sensitive to the thermal conductivity k of the crust and strongly suggest 
that K is large. We perform molecular dynamics simulations of the structure of the crust of an 
accreting neutron star using a complex composition that includes many impurities. The composition 
comes from simulations of rapid proton capture nucleosynthesys followed by electron captures. We 
find that the thermal conductivity is reduced by impurity scattering. In addition, we find phase 
separation. Some impurities with low atomic number Z are concentrated in a subregion of the 
simulation volume. For our composition, the solid crust must separate into regions of different 
compositions. This could lead to an asymmetric star with a quadrupole deformation. Observations 
of crust cooling can constrain impurity concentrations. 

PACS numbers: 97.60.Jd, 26.60.+C, 97.80.Jp, 26.50.+X 



I. INTRODUCTION 

What is the thermal conductivity of the crust of a neu- 
tron star? Recently the cooling of two neutron stars has 
been observed after extended outbursts [11, 0| ■ These out- 
bursts heat the stars' crusts out of equihbrium and then 
the cooling time is measured as the crusts return to equi- 
librium. The surface temperature of the neutron star in 
KS 1731-260 decreased with an exponential time scale 
of 325 ± 100 days while MXB 1659-29 has a time scale 
of 505 ± 59 days 0. Comparing these observations, of 
rapid coohng, to calculations by Rutledge et al. Q and 
Shternin et al. [3] strongly suggest that the crust has a 
high thermal conductivity. This would be expected if the 
crust is a regular crystal. 

In contrast, a low crust thermal conductivity, that 
would be expected if the crust is an amorphous solid, 
could help explain superburst ignition. Superbursts are 
very energetic X-ray bursts from accreting neutron stars 
that are thought to involve the unstable thermonuclear 
burning of carbon ^, Jl]. However, some simulations do 
not reproduce the conditions needed for carbon ignition 
because they have too low temperatures [7|. A low ther- 
mal conductivity could better insulate the outer crust 
and allow higher carbon ignition temperatures. 

The thermal conductivity is dominated by heat con- 
duction by electrons and this is limited by electron-ion 
scattering Q . Therefore in this paper, we present molec- 
ular dynamics simulations of the crust in order to calcu- 
late electron-ion scattering. We include many impurities 
based on results of a rapid proton capture nucleosynthesis 
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simulation [9| followed by calculations of electron capture 
[lol |. We find a regular crystal structure. We do not find 
an amorphous phase. We calculate the static structure 
factor S{q), that describes electron-ion scattering, and 
from S{q) we determine the thermal conductivity. Impu- 
rities can limit the thermal conductivity. If the impuri- 
ties are weakly correlated than their effect on the thermal 
conductivity can be described by an impurity parameter 

Q [UJ, 

Q = (AZ)2= (Z2>-(Z>2. (1) 

This depends on the dispersion in the char ge Z of each 
ion. The rp process ash composition of ref. [10| and ref. 
[l^ has a relatively large value oi Q = 38.9. Impurity 
scattering can be important at low temperatures where 
there is small scattering from thermal fluctuations. Note 
that ref. pH assumes the impurities are weakly corre- 
lated. If there are important correlations among the im- 
purities, for example if there is a tendency for low Z ions 
to cluster together instead of being distributed at random 
throughout the lattice, then the effects of impurities on 
the thermal conductivity could be different from what is 
calculated in ref. (llj . In this paper we perform MD sim- 
ulations to study the distribution of impurities and their 
effect on the conductivity. 

If the thermal conductivity is high, one may need ad- 
ditional heat sources in the crust in order to explain su- 
perburst ignition. Although Gupta et al. find some 
heating from electron captures to excited nuclear states, 
simple nuclear structure properties may provide a natu- 
ral limit to the total heating from electron captures Ts^ . 
Horowitz et al. fTi*] find additional heating from fusion 
of neutron rich light nuclei such as ^'^O+^'^O at densities 
near 10^^ g/cm"^. These fusion reactions are an important 
area for future work. Alternatively chemical separation 
with freezing, that was found in ref. could enrich 

the neutron star ocean with low Z elements and make it 
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enriched in low Z elements compared to the newly formed 
TABLE I : Abundance (by number) of chemical element Z. ,.j wix, j. j u • i i.- c ^^ 

^ ^ ' solid. What does chemical separation mean tor the struc- 

^ t^'^^ ^■'^^ crust? Perhaps the most conservative possi- 

bility is the following steady state scenario. We assume 
material accretes at a constant rate. Initially, chemical 

14 0023 separation enriches the ocean in low Z elements. Eventu- 

15 0023 ^lly ocean becomes so enriched that the composition 
20 0.0046 of low Z material in the newly forming solid is equal 
22 0.0810 to that in the accreting material. The system reaches 
24 0.0718 a steady state. The rate of low Z material accreting 

26 0.1019 into the ocean is equal to the rate freezing out (modulo 

27 0.0023 nuclear reactions). The sole effect of chemical separa- 
tion is to greatly enrich the ocean in low Z material. 
If we assume steady state, the composition of the crust 
will be the same as that of the original accreting mate- 

34 3866 Therefore, in this paper we perform MD simulations 

36 0023 determine the structure and thermal conductivity of 

47 0.0023 crust with the original Gupta et al. rp ash composition. 
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easier for superburst ignitition. 

In section In] we describe our molecular dynamics sim- 
ulations and the calculation of the thermal conductivity. 
Results for the structure of the crust, the static structure 
factor 5(9), and the thermal conductivity are presented 
in section ITTll We conclude in section Hvl 



II. MOLECULAR DYNAMICS SIMULATIONS 

In this section we describe our molecular dynamics sim- 
ulations and how we calculate the thermal conductivity. 
We begin with a discussion of our initial composition. 

A. Compositon 

Our model for the composition of the crust is the same 
as was used in previous work on chemical separation 
when the crust freezes [13] . Schatz et al. have calcu- 
lated the rapid proton capture (rp) process of hydrogen 
burning on the surface of an accreting neutron star [9|, 
see also This produces a variety of nuclei up to 

mass A « 100. Gupta et al. [I^ then calculate how 
the composition of this rp process ash evolves, because 
of electron capture and light particle reactions, as the 
material is buried by further accretion. Their final com- 
position, at a density of 2.16 x 10^^ g/cm"^, has forty % of 
the ions with atomic number Z ~ 34, while an additional 
10% have Z = 33. The remaining 50% have a range of 
lower Z from 8 to 32. In particular about 3% is ^^O and 
1% ^^Ne. This Gupta et al. composition |10|] is listed in 
Table HI In general, nuclei at this depth in the crust are 
expected to be neutron rich because of electron capture. 

Material accretes into a liquid ocean. As the den- 
sity increases near the bottom of the ocean, the material 
freezes. However we found chemical separation when the 
complex rp ash mixture freezes (l^ . The ocean is greatly 



B. Simulations of Crust Structure 

In order to calculate the thermal conductivity of a mul- 
ticomponent system one needs to understand its state. 
Monte Carlo simulations [lB| of the freezing of a classical 
one component plasma (OCP) indicate that it can freeze 
into imperfect body centered cubic (bcc) or face-centered 
cubic (fee) microcrystals. Unfortunetly not much has 
been published on the freezing of a multi-component 
plasma (MCP). There are many possibilities for the state 
of a cold MCP TJj. It can be a regular MCP lattice; or 
microcrystals; or an amorphous, uniformly mixed struc- 
ture; or a lattice of one phase with random admixture of 
other ions; or even an ensemble of phase separated do- 
mains. We perform classical MD simulations to explore 
the state of our MCP solid. 

The electrons form a very degenerate relativistic elec- 
tron gas that slightly screens the interaction between 
ions. We assume the potential Vij{r) between the ith 
and jth ion is. 



[r) 



(2) 



where r is the distance between ions and the electron 
screening length is Ae = 7r^/^/[2e(37r^ne)^/'^]. Here is 
the electron density. Note that we do not expect our re- 
sults to be very sensitive to the electron screening length. 
For example, the OCP melting point that we found in 
ref. , using a finite Ae , agrees well with the result for 
Ae = 00. 

To characterize our simulations , we define an average 
Coulomb coupling parameter F for the MCP, 



F = 



(Z5/3)(^)l/3, 



(3) 



where the mean ion sphere radius is a = 

(3/47rn)i/3 a,nd 

n = rie/ {Z) is the ion density. The OCP freezes at F = 
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175. In ref. [l^l we found that the impurities in our MCP 
fowered the melting temperature until F = 247. Finally, 
we can measure time in our simulation in units of one 
over an average plasma frequency Wp, 



LOr, 



E 



Mi 



(4) 



where Mj is the average mass of ions with charge Zj 
and abundance Xj (by number). Note that there will 
be quantum corrections to our classical simulations for 
temperatures significantly below the plasma frequency. 



C. Thermal conductivity 

The thermal conductivity k has been discussed by 
Potekhin et al. [3]. We assume k is dominated by heat 
carried by electrons [H], 



T^klTue 
3m!; v 



(5) 



where the effective electron mass is m* = ep = [kp + 
TTie)^/^ with kp the electron Fermi momentum and nie 
the electron mass. The electron collision frequency v is 
assumed to be dominated by electron-ion collisions [sl], 
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Here a is the fine structure constant and A is the 
Coulomb logarithm that describes electron-ion collisions 
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Here e{q, 0) is the dielectric function due to degenerate 
relativistic electrons, [3], see Eq. 2.3 of [19]. Note that 
for simplicity we neglect second and higher Born correc- 
tions to electron ion scattering in Eq. [71 see for exam- 
ple [l^. We are interested in the difference in thermal 
conductivity for different solid structures. Second and 
higher Born corrections should be the same for the dif- 
ferent structures. Finally, the lower limit qo in Eq. [7] is 
in a liquid phase and go = (67r^n)^/^ in a crystal phase 
i- 

The static structure factor S{q) describes electron- 
ion scattering. We calculate S{q) directly as a density- 
density correlation function using trajectories from our 
MD simulations, 



5(q) = (p*(q)p(q))-|(p(q))P 
Here the charge density p{q) is. 



(8) 



(9) 



with N the number of ions in the simulation and Zi, 
are the charge and location of the ith ion. We evaluate 
the thermal average in Eq. [8] as a time average during 
our MD simulations. 

The static structure factor S{q) includes both Bragg 
scattering contributions from the whole crystal lattice 
S'bragg (<z) and inelastic excitation contributions S'{q) 0], 



5(q) - S'i<D + 5bragg(q) 



(10) 



The Bragg contribution is a series of delta functions at 
momenta related to one over the lattice spacing. This 
describes Bragg scattering and helps determine the elec- 
tron band structure. It does not limit the electron mean 
free path. Instead the mean free path and thermal con- 
ductivity are determined by 5'(q). 

Our MD simulations are classical. Unfortunately this 
classical approximation makes the separation of S{q) into 
S' and "S'bragg somewhat ambiguous. We approximate 
5'(q) with a simple numerical filter applied to S{q). The 
filter removes delta function like contributions to S{q) 
that have a very rapid q dependence, and also removes 
numerical noise. This is discussed further in Section [IIII 



III. RESULTS 

To explore possible states for the multicomponent 
plasma we perform two molecular dynamics simulations. 
The initial conditions of these simulations are similar to 
those in [l^. The composition is indicated in TableHl We 
start by freezing a very small system of 432 ions. Here 
the ions were started with random initial conditions at a 
high temperature T and T was reduced in stages (by re- 
scaling velocities) until the system freezes. For the first 
simulation run, called rpcrust-Olb in Table [TTl we place 
four copies of this 432 ion solid in a larger simulation 
volume along with four copies of a 432 ion liquid con- 
figuration. This 3456 ion configuration is evolved at a 
lower temperature until the whole system freezes. Next, 
we evolve the 3456 ion solid at a reference density of 
n = 7.18 X 10~^ fm~^ (or 1 x 10^^ g/cm^) for a total 
simulation time of 2.4 x 10^ fm/c (8.9 x 10^ lj-^). The 
temperature was started at 0.325 MeV and slowly de- 
creased to a small value by the end of this time. The 
density and initial temperature correspond to F = 261.6. 
Evolution was done using the velocity verlet algorithm 
[20| using a time step of At = 25 fm/c for a total of 
9.6 X 10^ steps. This took about 2 months on a single 
special purpose MDGRAPE-2 [U board. Next, this low 
temperature configuration was reheated to T = 0.325 
MeV and evolved for 1.6 x 10^ fm/c. The total time was 
4 X 10^ fm/c. This somewhat complicated procedure was 
done for historical reasons. It does allow plenty of time 
for ions to diffuse throughout the simulation volume. 

Note that at our artificially high reference density (10^^ 
g/cm'^) free neutrons will be present. However, we are 
primarily interested in lower densities with out free neu- 
trons. Our results can be scaled to other densities and 
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TABLE II: Computer Simulations. The start time is ti, the 
finish time tf , N is the number of ions, and the temperature 
is T. Each simulation is at a density n = 7.18 x 10~^ fm""^ 
(1 X 10^^ g/cm^). Note that one over the plasma frequency is 
1/iOp = 270 fm/c. 

Run iV ti(fm/c) t/(fm/c) r(MeV) 

rpcrust-Olb 3456 1.6 x 10^ 0.325 

rpcrust-05 3456 4 x 10* 0.1 

4 x 10* 8 X 10* 0.2 

8 X 10* 1.2 X 10^ 0.3 

OCP 1024 1.6 X lO'^ 0.334 



temperatures such that the Coulomb parameter F re- 
mains the same, see below. Furthermore, although we 
quote all simulation times in fm/c, the times can be ex- 
pressed in terms of one over the average plasma frequency 
using 1/ujp = 270 fm/c. 

The initial configuration for run rpcrust-Olb, see Table 
im is shown in Fig. [T] The system is seen to be composed 
of two micro-crystals of different orientations. This is 
similar to the micro-crystals found in ref. [l^ upon freez- 
ing a one component plasma. In Fig. [T]we highlight the 
positions of the ^"'O ions (as small red spheres). These 
ions are located both in the crystal planes and in between 
them. The O ions are not spread uniformly throughout 
the volume but there is a tendency for them to cluster. 
This will be discussed in more detail below. 

This configuration was then reheated to T = 0.325 
MeV and evolved for 1.6 x 10^ fm/c. The final config- 
uration of run rpcrust-Olb is shown in Fig. (5] The two 
micro-crystals of different orientation are now gone. The 
system has managed to anneal into a single crystal with 
a single orientation. This suggests that micro-crystals 
could be an artifact of computer simulations of limited 
size and duration. It also suggests that neutron star crust 
could be formed with relatively large domain sizes. 

Figure [2] shows that O ions and other low Z impuri- 
ties are enhanced in regions on the left and right of the 
simulation volume. Because of the periodic boundary 
conditions this actually corresponds to a single region. 
We conclude that this complex mixture does not form a 
single uniform solid phase. Instead it separates into two 
solid phases. One phase is enriched in high Z ions and 
the other phase is enriched in low Z ions. 

To study this further we have performed another sim- 
ulation labeled rpcrust-05 in Table HIl The starting point 
was similar to run rpcrust-Olb with eight copies of a 432 
ion configuration placed into a larger simulation volume. 
This 3456 ion configuration was evolved for 2.5 x 10^ 
fm/c as the temperature was slowly decreased from 0.35 
MeV to a small value. Next, this low temperature con- 
figuration was heated to T = 0.1 MeV and evolved for 
400 million fm/c, the system was then heated to T ~ 0.2 
MeV and evolved for another 400 million fm/c and finally 
the system was heated to T = 0.3 MeV and evolved for a 
final 400 million fm/c as indicated in Table HIl The total 



simulation time including both the original preparation 
and the T = 0.1, 0.2, and 0.3 MeV runs was 3.7 x 10^ 
fm/c. 

The final configuration of run rpcrust-5 is shown in Fig. 
[3l The system involves only a single body-centered cubic 
(bcc) crystal. However O and other low Z ions are not 
uniformly distributed. Instead they are strongly enriched 
in a local region. This is indicated in Fig. 2] that shows 
the radial distribution function g{r) for run rpcrust-05 
at a temperature T — 0.1 MeV. The peaks in the Se- 
Se correlation function show the regular lattice planes. 
However g{r) for 0-0 is seen to be larger than one over 
a range of moderate distances r. This shows that the O 
ions are concentrated in a localized sub- volume. We con- 
clude that the complex rp ash mixture does not form a 
single solid phase. Instead, for this composition, the neu- 
tron star crust must be composed of two or more regions 
of different compositions. This disproves our steady state 
assumption. There appears to be no composition of the 
liquid ocean, no matter how enriched in low Z ions, that 
allows a uniform solid phase to form. 

These multiple regions of the crust with different com- 
positions may be very important for the structure of the 
neutron star. For example, if the phases are not dis- 
tributed uniformly, this could lead to a mass quadru- 
ple moment that might radiate gravitational waves psj . 
This nonuniform distribution of phases could arise from 
an anisotropic temperature because phase separation is 
temperature dependent. 

Finally for comparison we have also performed a one 
component plasma simulation, see run OCP in Table HIl 
where each ion has a charge Z = 29.4 equal to the av- 
erage charge in the MCP simulations. Simulation OCP 
started from a random configuration of 1024 ions and the 
temperature was reduced in stages until F = 300 at which 
point the simulation was observed to freeze. Finally this 
solid was heated up to F = 250 for the final results. 




FIG. 1: (Color on line) Configuration of the 3456 ion mixture 
in run rpcrust-Olb at the start of the simulation. The small 
red spheres show the positions of ^"^O ions, while ions of above 
average Z are shown as larger blue spheres. Finally, ions of 
below average Z (except for O) are shown as small white 
spheres. The left and right halves of the figure show two 
micro-crystals of different orientations. 
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FIG. 2: (Color on line) Configuration of the 3456 ion mixture 
in run rpcrust-Olb after a simulation time of 1.6 x 10^ fm/c. 
The small red spheres show the positions of ions, while 
ions of above average Z are shown as larger blue spheres. 
Finally, ions of below average Z (except for O) are shown as 
small white spheres. 
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FIG. 4: (Color on line) Radial distribution function g(r) ver- 
sus r over the mean ion sphere radius a, for run rpcrust-05 at 
a temperature T = 0.1 MeV. The dashed (red) line shows the 
correlation function between O ions while the solid (black) 
line shows the Se-Se correlation function. 



TABLE III: Coulomb Logarithm, Eq. [T] for a one component 
plasma (OCP). The Aocpflt value is from the S' {q) fit in ref. 
[31 while Aocp is our calculation for simulation OCP. 

r Apcpflt Aqcp 

250 0.362 0.348 



FIG. 3: (Color on line) Configuration of the 3456 ion mixture 
in run rpcrust-05 after a simulation time of 1.2 x 10^ fm/c. 
The small red spheres show the positions of ions, while 
ions of above average Z are shown as larger blue spheres. 
Finally, ions of below average Z (except for O) are shown 
as small white spheres. The concentration is seen to be 
enhanced in a sub-region to the right of center. 



A. Static Structure Factor 

We calculate the static structure factor S'(g) from the 
density-density correlation function, Eq. [51 The ther- 
modynamic average is approximated as a time average 
over 6.25 x 10^ fm/c of simulation time. We present 
results for the angle averaged S(ci) after averaging over 
approximately 50 different directions of q. These re- 
sults are somewhat time consuming because we calculate 
S{q) for approximately 1400 different values of \q\ for run 
rpcrust05. 

We calculate the inelastic contribution S'{q) by ap- 
plying a simple numerical filter that removes very rapid 
changes in 5'(g) with q. Our filter, applied to a table 
of qi and S'(qi) values, works as follows: if Sijqi) differs 
by more than some threshold ~ 0.1 from S{qi-i) than 



qi and S{qi) are removed from the table. This removes 
numerical noise and may remove delta function like con- 
tributions from the Bragg peaks. In addition, we may 
simply miss some Bragg peaks because we only calculate 
S{q) for a finite number of q points. Our motivation for 
this simple procedure is to calculate S'{q) and A based 
on S{q) calculations that are not likely contaminated by 
Bragg contributions. 

We first test this procedure with the one component 
plasma simulation OCP of Table [Til see Fig. [51 Our re- 
sults for S'{q) show more structure than the simple fit 
presented in ref. [H . Note that this may reflect a limita- 
tion of the fit. In addition, there is some high frequency 
noise in our simulation. However, there is good agree- 
ment, to 4 %, between our OCP simulation and the fit 
for the integral of S'{q) over q that is needed to calcu- 
late the Coulomb logarithm A, see Eq. [3 and Table IIIII 
Therefore, our procedure for S'{q) reproduces the known 
Coulomb logarithm and thermal conductivity of a one 
component plasma. 

Figuresini [3 and[51show S{q) for run rpcrust-05 at tem- 
peratures of T = 0.1, 0.2, and 0.3 MeV respectively. We 
expect similar results for run rpcrust-Olb. These figures 
also show the simple fit to S'{q) for an OCP presented 
in ref. [s!]. This fit is significantly below S'{q) for run 
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FIG. 5: (Color on line) Inelastic contributions to the static 
structure factor S'{q) versus momentum transfer q times the 
mean ion sphere radius a, for run OCP, solid black line and 
the simple fit presented in ref. [8j , dashed red line. 



TABLE IV: Coulomb Logarithm, Eq. [T] for rp process ash 
composition. The A values are from run rpcrustOS at the 
indicated temperatures T, while Aocpflt is from ref. for 
a pure OCP and Aocp+imp also includes impurity scattering 
from ref. with Q = 38.9. 

T(MeV) r Aocpflt Apcp+imp A 
0.1 850 0.104 0.146 0.173 
0.2 425 0.232 0.276 0.366 
0.3 283 0.334 0.377 0.530 



rpcrust-05. Finally, these figures show the contribution 
of impurity scattering from [llj added to the OCP fit re- 
sults. Impurity scattering depends on Q, see Eq. [T]and 
Q = 38.9 for run rpcrust-05. We find that results for run 
rpcrust-05 are still above the OCP fit even when impurity 
scattering is added. Note, that impurity scattering is au- 
tomatically included in our MD simulation because of the 
complex composition used. Table IIVI presents Coulomb 
logarithms A for rp ash composition. Again, results for 
run rpcrust-05 are above the OCP plus impurities calcu- 
lation. However the difference is only 18% at a tempera- 
ture of 0.1 MeV. Note in Reference yuj it was explicitly 
assumed that the impurities are randomly distributed. 
However, we find strong correlations among the impuri- 
ties, see Fig. 2] for example. Therefore it is perhaps not 
surprising that we find larger effects from impurities than 
ref. [ni. 

Our results for S{q) in Figs. [5][8]show statistical noise. 
However some of the effects of this noise average to zero 
when one integrates over S'{q) to calculate A. We es- 
timate the statistical error in our calculation of A at 



FIG. 6; (Color on line) Static structure factor ^(g) for run 
rpcrust-05 versus momentum transfer q times the mean ion 
sphere radius a, dotted black line, at a temperature T = 0.1 
MeV. The solid black line is an approximation to the inelastic 
contribution S'{q). This is calculated with a simple numerical 
filter applied to S{q). Finally the dashed green line is the fit 
to OCP results for S'{q) from ref. and the dashed dotted 
red line adds impurity scattering from ref. [ll|] to these OCP 
results. 



TABLE V: Thermal conductivity for a temperature of T = 
0.043 MeV (5 x lO^K). Results have been scaled to the in- 
dicated densities. The thermal conductivity k is from run 
rpcrust-05 while Kocp+imp is based on the OCP plus impu- 
rity values Aocp+imp in Table IIVI and kqcp is for the fit to 
an OCP without any impurity scattering. 

r p KOCP KOCP+imp K 

(g/cm^) (erg/K cm s) (erg/K cm s) (erg/K cm s) 

850 7.91 X 10" 2.48 x lO^'' 1.77 x 10^" 1.49 x 10^^ 
425 9.89 X 10"' 5.58 x 10^® 4.69 x 10^** 3.54 x 10^* 
283 2.92 X 10^" 2.58 x 10^* 2.29 x 10^** 1.63 x 10^* 



T = 0.1 MeV, see Table JVl to be ± 0.001 by compar- 
ing calculations of A using configurations for simulation 
times of 3x 10* fm/c to 4x 10® fm/c to a calculation using 
configurations from 2x 10* to 3x 10* fm/c. We emphasize 
that our procedure to calculate S'{q) from S{q) is model 
dependent. Our numerical filter not only removes Bragg 
peaks but it may also remove some statistical noise. Note 
that removing some noise seems to have minimal effects 
on the values of A that we calculate. We do not believe 
our results in Table HVl are very sensitive to our procedure 
to determine S'{q). This is based on explicit calculations 
with a few different procedures. 
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FIG. 7: (Color on line) Static structure factor S{q) for run 
rpcrust-05 versus momentum transfer q times the mean ion 
sphere radius a, dotted black line, at a temperature T = 0.2 
MeV. The solid black line is an approximation to the inelastic 
contribution S'{q). This is calculated with a simple numerical 
filter applied to S{q). Finally the dashed green line is the fit 
to OCP results for S'{q) from ref. and the dashed dotted 
red line adds impurity scattering from ref. [ill ] to these OCP 
results. 



B. Thermal Conductivity 

We now calculate the thermal conductivity k using our 
results for the Coulomb logarithms. These results can be 
scaled to a range of densities n and temperatures T so 
that the value of F, Eq. [31 remains the same. Table fVl 
presents k at a temperature of T = 5 x 10® K (a typical 
value for a super bursting star). The thermal conduc- 
tivity is lower for run rpcrustOS than for an OCP. First, 
this is because run rpcrust05 has a large number of im- 
purities, corresponding to the large impurity parameter 
Q = 38.9. Second, we think k may be further reduced be- 
cause the impurities in run rpcrustOS are not distributed 
uniformly. Instead they are concentrated in one region. 

Although our simulations show some of the effects of 
impurities on the thermal conductivity, we emphasize 
that there may be important finite size effects because 
we find clustering. It is unreahstic to describe a large 
system by simply repeating our small simulation volume 
many times. This would describe the impurities as being 
concentrated into many very small regions. Instead, we 
believe the concentration of impurities indicates phase 
separation. We think that a large sample will separate 
into two (or more) bulk phases. It is important to study 
phase separation further with larger molecular dynamics 
simulations and this may change our thermal conductiv- 
ity results. In general, one phase will be enriched in high 
Z ions while the other is enriched in low Z ions. Phase 
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FIG. 8: (Color on line) Static structure factor S{q) versus 
momentum transfer q times the mean ion sphere radius a for 
run rpcrust-05, dotted black line, at a temperature T — 0.3 
MeV. The solid black line is an approximation to the inelastic 
contribution S'{q). This is calculated with a simple numerical 
filter applied to S{q). Finally the dashed green line is the fit 
to OCP results for S'{q) from ref. and the dashed dotted 
red line adds impurity scattering from ref. [ll|] to these OCP 
results. 



separation may act to reduce the impurity parameter Q 
and increase the thermal conductivity. For example, Q 
will be reduced in the high Z phase because low Z im- 
purities have gone into the other phase. 

In addition, nuclear reactions may reduce Q further. 
In general, we expect nuclear reactions to preferentially 
burn low Z impurities because of their low Coulomb bar- 
riers. See for example ref. [l3|- This will reduce Q and 
increase the thermal conductivity. One should study how 
Q evolves with depth because of reactions. Finally, it is 
important to analyze observations of crust cooling after 
extended outbursts [1, 0] to see what observational con- 
straints can be placed on the thermal conductivity and 
Q. It may be that observations of rapid crust cooling can 
strongly limit the size of Q [221 ■ 

Phase separation may have another important effect. 
It will create layers in the crust of different compositions 
and densities. These layers may not be spherically sym- 
metric. For example, phase separation depends on tem- 
perature. Therefore an anisotropic temperature distribu- 
tion will lead to an anisotropic density. It is important 
to study how phase separation will change the structure 
of the star. 
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IV. SUMMARY AND CONCLUSIONS 

The crust of an accreting neutron star, likely, has a 
complex composition with many impurities. Nuclei are 
synthesized via the rapid proton capture process and the 
composition is modified by electron capture as material is 
buried to greater densities. We have performed MD sim- 
ulations, with a complex composition, to study the struc- 
ture of the crust. Our simulations form ordered crystals 
rather than an amorphous solid. 

However, we find phase separation. Some low Z impu- 
rities are concentrated into a subregion of the full simu- 
lation volume. This phase separation, between two solid 
phases, is similar to the chemical separation found pre- 
viously between liquid and solid phases [12] . Previously, 
we assumed a steady state equilibrium where chemical 
separation greatly increases the concentration of low Z 
impurities in the liquid ocean. However, the composition 
of the solid crust was assumed to be the same as that 
of the accreting material. Our new results disprove this 
steady state assumption. 

The crust can not be uniform, given our initial com- 
position. Phase separation will divide the crust into two 
or more regions of different compositions. This may have 
important implications for the structure of the star. For 
example, composition anisotropies could lead to gravi- 
tational wave radiation from a quadrupole deformation 



[2 3 *1 . In future work we will study the size of possible com- 
positional asymmetries because of an anisotropic temper- 
ature distribution 

We calculated the static structure factor S{q) for our 
simulations and from the thermal conductivity k. 
Since our simulations have a complex composition, we 
automatically include the contributions of impurity scat- 
tering. We find that k is somewhat reduced because of 
impurity scattering and because the impurities are not 
distributed uniformly. We expect the same results for 
the electrical conductivity a, that is important for mag- 
netic field decay [26), and the shear viscosity ?7, that can 
damp neutron star oscillations [11], [2^. The reduction 
in K may be observable in crust cooling times and these 
observations may set limits on impurities. Future work 
should study how phase separation and or nuclear reac- 
tions impact impurity concentrations. 
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